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We consider the dispersion properties of tracer particles moving in non-equilibrium heterogeneous 
periodic media. The tracer motion is described by a Fokker-Planck equation with arbitrary spatially 
periodic (but constant in time) local diffusion tensors and drift, eventually with the presence of 
obstacles. We derive a Kubo-like formula for the time dependent effective diffusion tensor valid 
in any dimension. From this general formula, we derive expressions for the late time effective 
diffusion tensor and drift in these systems. In addition, we find an explicit formula for the late 
finite time corrections to these transport coefficients. In one dimension, we give closed analytical 
formula for the transport coefficients. The formulas derived here are very general and provide a 


straightforward method to compute the dispersion 
advection-diffusion systems. 
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I. INTRODUCTION 

Fokker-Planck (FP) equations, and their associated 
stochastic differential equations (SDE) arise in a huge 
variety of physical systems m- In fluid mechanics and 
hydrology, such equations describe the transport of pas¬ 
sive scalars advected by a fluid and locally dispersed by a 
microscopic molecular diffusivity. They also describe the 
diffusion of colloids and polymers in soft matter physics, 
as well as charge and energy transport in solid state 
physics. In many systems, damping due to the environ¬ 
ment is very strong and inertial effects are thus negligible. 

In such systems, the evolution of the probability density 
function (pdf) p(x, t) of a tracer particle at position x 
and time t obeys the FP equation 

dtp = dxi [dxj [Kij{x)p] - Mi(x)p} = -i/xP, (1) 

where the Einstein convention on indexes summation is 
used, and the indexes i,j run from 1 to the spatial di¬ 
mension d. The time independent tensor Kij (x) is the lo¬ 
cal diffusivity tensor and can, in general, vary in space in 
media of non-uniform composition, since diffusivity is de¬ 
pends on the local material properties of the system. The 
local diffusivity may also vary in space due to tempera¬ 
ture inhomogeneities, or in the case of particle diffusion 
in viscous fluid flows close to solid boundaries [1]. The 
term Ui(x) represents a drift term generated by forces 
such as buoyancy forces or electromagnetic forces acting 
on colloids, fluid flows, as well as thermodynamic effects, 
such as thermodiffusion in the presence of temperature 
gradients . The presence of impenetrable obstacles can 
also be treated by introducing surfaces at which there is 
no normal flux, 

ni{x){dx,j[Kij{K)p]-Ui{x)p} = 0, xeS, (2) 

where rii represents the normal vector to the surface S 
of the obstacles. 

Equation Q describes the motion of tracer particles at 
small scales in a wide range of heterogeneous media. At a 


properties in arbitrary non-equilibrium periodic 


coarse-grained level, the motion is described by effective 
transport coefficients characterizing the average flux of 
particles and the relative spreading with time of initially 
close particles. From a theoretical point of view, the 
dispersion of tracer particles is defined in terms of two 
quantities, the mean displacement Ri(t) during a time t 
in the direction i 

/r,(t) = (X,(t)-X,(0)), (3) 

(where X(t) denotes the position of a tracer particle and 
(• • •) denotes ensemble averaging), and the correlation of 
the dispersion in the directions i and j defined by 

a.,(t) = (4) 

The quantity crij{t) thus characterizes the dispersion of 
a cloud of particles about its mean position. Of partic¬ 
ular interest is the large time behavior of the dispersion 
properties, 

/ii(t) ~ 14 t, (t-)■ oo) (5) 

~2 t {Dij F Cij/t), (t-)-oo), (6) 

where 14 is the average particle drift, Dtj is the late time 
effective diffusion tensor, and Cij gives the leading order 
corrections and provides information on how fast the sys¬ 
tem approaches the diffusive limit. We will see that the 
asymptotic form (§ is valid if the medium is spatially 
periodic. The quantities 14 and Dij characterize the late 
time, large scale transport properties of the system and 
are routinely measured in experimental systems by sin¬ 
gle particle tracking methods [5] or by ensemble based 
measurements such as fluorescence recovery after photo- 
bleaching (FRAP) [7]. Effective transport coefficient are 
also important for estimating the spread of pollutants, 
filtration, mixing and chemical reaction times [5]. 

Transport and dispersion properties in non-equilibrium 
periodic systems have been investigated at length in var¬ 
ious contexts. In a pioneering, classical work, Taylor [9] 
considered the dispersion of particles moving in viscous 
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flows in cylindrical containers, and showed that the ef¬ 
fective diffusion coefficient could be orders of magnitudes 
larger than the molecular diffusivity. Effective transport 
in incrompressible fluid flows were also investigated at 
length, such as in the case of diffusion in Rayleigh-Benard 
convection cells [TOHE], in frozen turbulent flows [13] or 
in porous media [14H20j using homogenization theory. 
However, these approaches have been restricted so far 
to incompressible flows and it is not clear how to gener¬ 
alize them. In another context, results were also derived, 
with the methods of statistical physics, for particles dif¬ 
fusing in periodic potentials when the diffusivity is con¬ 
stant pTl425j : in here the potential barrier slows down 
the diffusion as the particle tends to become trapped in 
local mining. The non-equilibrium problem of transport 
in one-dimensional potential in the presence of external 
constant force (still with uniform microscopic diffusivity) 
was also considered |26H29]; this system is predicted to 
exhibit a huge increase of the effective diffusion near a 
critical tilting force. Results have also been obtained in 
the more general case where the noise amplitude is a pe¬ 
riodic function of position |5{H - I5?| , notably in the context 
of the modeling of transport in periodic channels within 
the Fick-Jacobs approximation (where the full Fokker- 
Planck equation is approximated by an effective one di¬ 
mensional one). In many of the above cases, results only 
exist in one dimension, however their generalization to 
higher dimensions is implicit in the general results we 
will present here. 


In this paper, we present a general formalism to com¬ 
pile the time-dependent dispersion tensor g^, (t) (section 


nil, the late time diffusion tensor Dij (section |rV| and 


the corrections Cj, describing the approach to the diffu¬ 
sive limit (section]^. A summary of the results is given 
in section]^ for the reader who might not be interested 
in the technical aspects of the calculations. The results 
derived are valid for any system described by Eq. 0 for 
periodic Kij(x) and Mi(x); they encapsulate all known re¬ 
sults for this problem and also generalize them to more 
complex situations. We will explicitly compare them to 
the expressions of the literature in various cases (section 


VI), and we will also give an explicit fully analytical so¬ 


lution for the one dimensional problem (section VII). A 
brief derivation of our main results has already appeared 
in Ref. |33j , where they were used to analyze diffusion in 
a system with a spatially varying diffusivity in the pres¬ 
ence of a uniform applied force. The present paper ex¬ 
tends this recently developed approach, new features are 
the computation of the time-dependent dispersion prop¬ 
erties (instead of only the late time diffusive limit) and 
that the incorporation of the effect of impenetrable obsta¬ 
cles with no-flux boundary conditions. Furthermore, the 
derivation presented here is different from that appearing 
in Ref. [33], as it is not based on the SDE formulation 
used there. Finally, the present paper contains explicit 
expressions for the first temporal corrections Cij to the 
late time diffusion tensor. We will show that the 1/t 
decay of the late time correction [Eq. 0 ] occurs in all 


dimensions. This universality can be understood as, for 
periodic systems, the Fokker-Planck equation restricted 
to the basic cell 17 will in general have a gap and not 
a continuum (whose density would be dependent on the 
spatial dimension) of eigenvalues near zero. 


II. NOTATIONS AND SUMMARY OF RESULTS 


We assume that the space can be divided into individ¬ 
ual cells U, and we that the diffusivity tensor Kiji'x.) and 
the drift Wi(x) are periodic in space, and thus satisfy 

KyXx) = Ky(x-I-k), ?.ti(x) = Mi(x-f k), (7) 

for all vectors k joining the centers of the cells compos¬ 
ing the periodic structure. Denoting by X(t) the position 
of an individual particle at time t, we define the propa¬ 
gator p(x,t|y) of the stochastic process X(t) in infinite 
space; p(x, t|y) is thus the probability density to observe 
a particle at position x at time t starting from y at time 
0. Mathematically, p is the solution of the equation in 
infinite space 

[dt + H^)p{x,t\y) = 0 ; p(x, 0|y) = (5(x - y). (8) 

We also introduce the process X(t) as the position of the 
particle X ’’modulo U”, X is defined as the position of the 
particle in a reference cell Uq obtained after an integer 
number of translations of X along the vectors joining the 
centers of the cells composing the periodic structure. We 
introduce the propagator P(x,t|y) for the process mod¬ 
ulo U, P thus represents the probability density, starting 
at position y at time 0, to be at time t at position x 
modulo U. P is the solution of the FP equation Q with 
periodic boundary conditions on the boundaries of U and 
initial value P(x, 0|y) = (5(x — y). There is an obvious 
relation between these propagators, 

t\y) p(x + k, t|y), (9) 

k 

where k represents the vectors joining the center of the 
cell containing x to the centers of all the other cells of 
the periodic structure. The value of P at late times is 
denoted by Pq and is the steady-state distribution of par¬ 
ticles in one period (modulo translations along the lattice 
vectors). In contrast, the propagator in infinite space p 
does not tend to a steady state, but becomes Gaussian 
at long times, with a variance growing linearly as t, rep¬ 
resenting the spreading of the particle density. In what 
follows, we always assume that at t = 0 we are in the 
steady state, that is to say that the process X(t = 0) 
modulo U has the distribution Pq. 

Before entering the details of the calculations, let us 
give briefly our main results. First, we will show that the 
temporal Laplace transform aij{t)e~^^dt of 
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the dispersion tensor is given by 

^ J Kij{x)Po{x) 

dxdy[ui{x)u* (y) + Uj (x)u*(y)]P'(x, s|y)Po(y)- 

( 10 ) 


«2 , 
s jjn 


In this relatively compact formula, P'(x, s|y) represents 
the Laplace transform of P(x,t|y) — Po(x). We have 
also defined an effective velocity field u that takes into 
account the presence of obstacles (if present) 


u^{x) = Mi(x) - nj(x)«;ij(x)i5s(x), 


( 11 ) 


where 6$ is a surface delta function; to be clear with the 
notations we precise that 

/ (ixn^(x)Ky(x)(5s(x)V'(x) = / dSjKij{x.)'tjj{x) (12) 
da Jsn 

for any function tp, with the infinitesimal normal surface 
vector dSi and the unit normal vector Ui oriented to¬ 
wards the interior of the obstacles. Finally, we have also 
introduced the drift field u* which can be interpreted as 
the drift field for the particles upon time reversal, and is 
given by 


'(x) = Mi(x) - 2 


Po(x)' 


(13) 


where we have introduced the current in the steady state 

J^ 

J,(x) = Mj(x)Po(x) - a 2 ,jKjfc(x)Po(x)]. (14) 

The notation u*(x) denotes that we have added surface 
drift terms which prevent the particles from entering the 
obstacles. 


((x) = u* (x) - nj(x)Kij(x)5s(x). 


(15) 


Equation ( |10[ ) states clearly that the large scale disper¬ 
sion properties can be obtained from the properties of the 
propagator (with periodic boundary conditions), which is 
used to average the drift fields with and without time- 
reversal; this is the first result of the paper. 

The second main result is an expression for the late 
time diffusion tensor Dij. Taking the large time limit of 
Eq. (10) leads to 


Dij = dx Ky (x)Po(x) 

Jn 

dx [ui{x)fj{x) +Uj{x)f,{x)] 

dSi{x) [kj/(x)/,(x) + Kj 7 (x)/,(x)], (16) 

where the function fj is defined by 

/j(x) =- / dy Po(y)G(x|y)M*(y) (17) 
Jn 


1 

2 JO 

i f 


Sn 


and 

poo 

G'(xly) = / dt[P(x,t|y) - Po(x)]. (18) 

Jo 

G is the pseudo-Green function [M] of the operator H. 
For computational purposes, it is useful to note that fi 
can be obtained by solving the non-homogeneous linear 
partial differential equation 

= - u*(x)Po{x) + Po(x) [ dy Ui{y)Poiy) 

Jn 

- Po{x) [ dSi{y)Po{y)Ku{y) (19) 

JSn 

with periodic boundary conditions and also with the con¬ 
dition 

nj{uj{x)f,{x) - d^^[Kkj{x.)f^{x)]} = -njKj,{x)Po{x) 

( 20 ) 

at the surface of the obstacles, together with the orthog¬ 
onality condition 


[ dx fi{x) = 0. 
Jn 


( 21 ) 


The above expressions are exact, general and can be used 
to compute the dispersion properties in any advection- 
diffusion problem with periodic diffusivity tensor and 
drift fields, in any dimension, in the presence of obsta¬ 
cles with reflecting boundaries. We will show that our 
expression for Dij contains existing expressions for the 
effective diffusion tensor in porous media for particles 
transported by incrompressible fluids [Ill [35], diffusion in 
periodic potentials |36j . and diffusion in non-equilibrium 
one dimensional systems HH |32| . 

One may also derive find an alternative expression for 
Dij by defining a function f* as 

fHy)=- [ dx Po{y)G{x\y)u^{x), (22) 

Jn 

in terms of which the effective diffusion tensor reads 

Ai = dy K^j{y)Po{y) 

Jn 

+ \ K {y)fj (y) + u*iy)f* (y)] 


'Sn 


dSi{y)[Kuiy)f^y^) + Kji{y)fi{y)]- (23) 


This function f* can also be obtained as the solution of 
the partial differential equations 

HyfHy) = - My)Po{y) + A(y) f dx [u*(x)Po(x)] 

Jn 

- Po{y) f dSi{x)Po{x)Kii{x). 

JSn 


(17) 


(24) 
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FIG. 1: (a) Examples of stochastic trajectories of two tracer 
particles diffusing in a heterogeneous medium (here a 2D ar¬ 
ray of circular obstacles) under a force F (arrow) that biases 
the motion to the right, (b) Elements of the dispersion tensor 
for the system desc ribed i n (a), calculated from our theoret¬ 
ical formulas Eqs. ( |16|19[ ) (full and dashed lines) vs results 
of stochastic simulations (symbols). Parameters: disk radius 
R — 0.251/ with L the length of the side of the periodic cell. 
Results for Dij are expressed in units of the (constant) diffu¬ 
sion coefficient of the tracers outside the obstacles. 


Although numerous expressions exist in the literature for 
Dij for particular cases, the correction tensor Cij has to 
our knowledge only been studied for diffusion in a peri¬ 
odic potential m or with periodic diffusivity [35]) the 
above expressions enable its computation in a system¬ 
atic way and at a computational cost equivalent to that 
required to calculate Dij. 

Our results can be applied to study tracer particles 
diffusing in a 2D array of circular impenetrable obsta¬ 
cles (hard disks) [FigQa)] and subjected to an external 
force. For this system, the elements of the effective diffu¬ 
sion tensor match perfectly with the results of stochastic 
simulations (the algorithm is presented in appendix [A| . 


III. DERIVATION OF THE KUBO FORMULA 
FOR THE TIME DEPENDENT DISPERSION 

Here we describe the derivation of the Kubo formula 
for the dispersion tensor. We first briefly compute the 
average displacement in the direction i during a time t, 
defined in Eq. (|^, which is given by 


Fi{t) = [ dx [ dy Po{y)p{x,t\y){xi - yi). (29) 
Jn 

Note that the integral over the whole space denotes 
the integral over the volume not occupied by obstacles 
(where p does not vanish). Taking the time derivative of 
the above expression, and using Eq. Q, we obtain 

dtFz{t)= [ dx [ dyPo{y){-d^^[uj{x)p{x,t\y)] 

Jn 

+ da;^da;^^[Kkj{x)p{x,t\y)]}{xi -y,). (30) 


where we have denoted H* the transport operator for the 
time-reversed process, 

Hy^-) = 9yAK{y)i-)] - dy,dy^[Kij{y){-)]. (25) 

The function f* has periodic boundary conditions, and 
satisfies at the surface of the obstacles 


Integrating by parts over x and applying the divergence 
theorem yields 

dtPzit) = dx dyPo{y){uAx)p{x,t\y) 

Js.<i Jn 

-d^^[Kikix)p{x,t\y)]}, (31) 


nj{y){uj{y)f*{y) - dy^[Kkj{,y)fHy)]} = -njKji{y)Po{y), 

(26) 

and also satisfies the orthogonality condition 

/ dy/;(y) = 0. (27) 

Jn 

Thus, we have two formulas for the effective diffusion 
tensor, which can be straightforwardly evaluated numer¬ 
ically by solving a set of elliptic partial differential equa¬ 
tion on the domain H. 

Einally, we will show that the tensor Cij describing the 
approach to the diffusive limit can be obtained from 


where we have taken into account the reflecting boundary 
conditions at the obstacles, so that no surface integral has 
appeared. Decomposing the integral over x in integrals 
over each unit cell, using the periodicity property ([^ of 
Ui and Kij and the definition (§, we get 


dtPi{t) 



dyPo{y){u^{x)P{x,t\y) 


-d^^[K,k{x)P{x,t\y)]}. 


(32) 


We recall that, since Pq is the steady state particle dis¬ 
tribution, it satisfies 


^ /■ , /)(x)/;(x)-k/j-(x)/*(x) 

*'“2^ Po(x) 


dyPo{y)P{x,t\y) = ^o(x). 


(28) 


(33) 
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Using this property, and again applying the divergence 
theorem, (32) becomes 


[ dS'^(x)Kij(x)Po(x)} = U, 
■JSa 


(34) 


where Sq represents the surface of the obstacles inside the 
cell U, and the infinitesimal surface vector dSj is oriented 


towards the interior of the obstacles. Eq. (34) shows 


that, when the initial conditions are those of the steady 
state, the average drift Vi is constant in time. This is in 
agreement with the classic result of Stratonovich [5^ . 

Let us now derive a Kubo formula for the time depen¬ 
dent dispersion tensor cry [defined in Eq. (|^]. Define ipij 
to be 

aij(t) = (35) 

By definition, the expression for ijjij is 

j dyPo{y)p{:>c,t\y){xi-yi){xj-yj). 

Jn 

(36) 

We denote by (j){s) = dt(j){t)e~^* the temporal 

Laplace transform of any function (p. The Laplace trans¬ 
form of the FP equation ([^ is given by 

sp(x,s|y) = 6{x-y) 

+ doo, {(92.J/«*j(x)p(x,s|y)] - Ui(x)p(x,s|y)} . (37) 

Using the above equality, the expression for ipij becomes 

'0it(s) = - / dx dy Po{y){x^ - yi){xj - yj) 

■s ./B'i Jn 

dxk {dxi [Kfci(x)p(x, sjy)] - ■Ufe(x)p(x, sjy)} . (38) 
We remark that 

dxk[{^i ~ ~ Vj)] — dki{xj — yj) 

+ Skj{x,-y,), (39) 


so that the integration by parts over x in Eq. (38) leads 
to an expression of the form 


i’ijis) = Eij{s) + Eji{s) 


(40) 


with 


1 


Eijis) =- dx dy Po{y)ixj - yj) 
s Jn 

{u^{x)p{x, sjy) - dx^ [Kki{x)p{x, sjy)]} . (41) 

Integrating by parts the term containing Kik leads to 

Etj{s)=^ [ dx [ dy Po{y){xj-yj)u,{x)p{x,s\y) 

s Jn 

- dSi{x) / dyPo{y)Kuix){xj-yj)p{x,s\y) 

■s Js Jn 

+ - [ dx [ dyPo{y)p{x,s\y)Kij{x). (42) 

s JR'i Jn 


Using the periodicity property of k and the Laplace trans¬ 
form of the steady state property (33), we obtain 


/ dx Ui{x)Po{x) 

[ dx [ 

In 

Jr<‘ Jn 


' ^ [ dx Po{x)Kii{x) = s (43) 

Jn 

where represents the average over Pg- Decomposing 
the integral appearing in (42) over x on all the individual 


cells of the periodic structure, and changing of variable 
X—J-x-l-k, y—J-y-l-kin each of them (where k is the 
lattice vector such that x -|- k is inside the cell D), and 
summing over all lattice vectors k again, we see that we 
can exchange the integration domains between x and y, 
leading to 


Ei]{.s) = - [ dxui{x)Bj{x,s) 
s Jn 

1 


dSi{x)Kii{x)Bj{x,s) + 

■S Jsn « 


where we have defined 

Bj{x, s)= dy p{x, s\y)Poiy){xj - y^). 


(44) 


(45) 


Now comes the key point of our derivation. In order to 
calculate Bj, we first consider another probability distri¬ 
bution q defined as 


9(y,i|x) = 


p(x,t|y)Po(y) 

^^o(x) 


(46) 


Using Bayes’ theorem, we see that q{y,t\x) is the proba¬ 
bility density at the position y at time 0 given that x is 
the position at the later time t; q is thus the propagator 
of the tracer particles under time reversal. The evolution 
of q satisfies 

dtqiy,t\x)Po{x) = - a3,Jui(x)Po(x)g(y,t|x)] 

+ dxBx^[Kij{x)Po{x)q{y,t\x)]. (47) 

Expanding all the derivatives and using the fact that 
PPq = 0, we find that q satisfies 


dtq{y,t\x) =[u*{x)dx 
where 


Ktj{^)dxidx^]q{y,t\x), (48) 


<(x) = 2 


dx^[Kij{x)Po{x)] 

^o(x) 
>/o*(x) 


- Mi(x), 


= ^^(x) - 2 


^o(x) ’ 


(49) 


is the drift of the time reversed process and where we 
have introduced the current in the steady state Jt 

Ji{x) = Ui{x)Po{x) - 5a;jKjfc(x)Po(x)]. (50) 
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If we interpret (48) as a backward Fokker-Planck equa¬ 
tion [3], we deduce that q also satisfies the associated 
forward Fokker-Planck equation 

dtq{y,t\x) = 

- 5y,K(y)9(y>i|x)] +5y,5y^.[Ky(y)g(y,t|x)]. (51) 

Therefore q is the propagator of fictive particles, moving 
in an effective drift field u* instead of Ui. Let us say a 
little bit more about the properties of q. From its defi¬ 
nition (46), we immediately see that the initial condition 
for q is 


g(y,0|x) = 5(x-y). 


(52) 


Given that no particles can flow in or out of the obstacles, 
the time reversed process must also have no-flux bound¬ 
ary conditions at the surface of obstacles. That is to say 
that 

0 = ni{u*{y)q{y,t\x) - dy^ [K„ (y)g(y,t|x)]} 

^({<(y)^o(y) -5j/,Kj(y)Po(y)]}p(x,t|y) 


^o(x) 


- Kij(y)Po{y)dy^p{x,t\y) 


(53) 


and as a consequence we see that niKij{y)dy^p{x,t\y) = 
0, which recovers an established boundary condition for 
the no-flux boundary condition in terms of the starting 
coordinate y [3]. It is interesting to note that the steady 
state current of the time reversed process is given by — Ji, 
i.e. exactly the opposite of the current of the original 
process. In addition, we see that for currentless steady 
states, where Ji = 0, we have that u* = Ui and thus 
the original and time reversed processes are statistically 
identical in that p(x,t|y) = g(x,f|y). Interestingly, for 
advection of a particle with constant molecular diffusivity 
by an incompressible flow, one can easily show that Pq = 
I/|n| where |r2| is the volume of the unit cell and that 
Ji is non-zero and given by Ji = Ui/\Q.\, consequently 
u* = —Ui ,i.e. the time reversed process has the opposite 
flow field to the original process. 

Eq. (45) may be rewritten in terms of q as 


where we have used the no-flux condition (this time of 
the time reversed process) at the obstacles boundaries. 
Integrating by parts the term containing and then 
switching back to the propagator p instead of q and in¬ 
voking the periodicity property, we obtain 

By{x,s)=-- [ dy Po{y)P{x,s\y)u*{y) 
s Jn 

+ - [ dSi{y)Po{y)P{x,s\y)Kij{y). (57) 


Inserting the above expression into Eqs. (44),(36) gives 
a formal expression for the tensor tjjij. Using the fact 
that P(x, s|y) ~ Po(x)/s for s —0, one can check that, 
for large times, ipij ~ fj,i(t)fj,j(t). It is now useful to 
introduce the difference between the propagator and its 
stationary value 


P'(x, t\y) = P(x, t\y) - Po(x), 
and thus has Laplace transform 

P'(x,s|y) = P(x,s|y) - yPo(x). 

In terms of P', we 


(58) 


(59) 


obtain from Eqs. (|3|),(|3^,(|^, 
and (57) the following final expression for the temporal 
Laplace transform of the dispersion tensor try (t) 

dij{s) = ^ [ dx Kij{x)Po{x) 
s Jn 

+ ^ I dx[ui{x)Kj{x, s) + Uj{x)Ki{x, s)] 
s Jn 

-\[ dSi{x)[Ku{x)Kj{x,s) + Kji{x)Ki{x,s)], (60) 
S JSa 


with 


Kj{x,s) = - f dy Po(y)P'(x,s|y)M*(y) 

Jn 

+ [ dSi{y)Po{y)P'{x,s\y)Kij{y). (61) 
JSn 

and we recall that u* is the effective drift field after 


Bj{x,s) = [ dy q{y,s\x)Po{x){xj -pj). (54) time reversal symmetry defined by Eq. (@. Equations 

Jw^ (l4^.((60l) and (loTl) eive a closed form exoression for the 


Laplace transforming the Fokker-Planck equation (51) 
for q yields 

S 9 (y,s|x) = (5(x-y) 

-dyAu*iy)qiy,s\x)]+dy^dy^[K^j{y)q{y,s\x)]. (55) 


Inserting the above equality into Eq. (54) and integrating 
by parts over y, we obtain 

Bj{x,s) = - [ dy Po(x){-M*(y)g(y,s|x) 

+ dy^[Kkj{y)q{y,s\x)]}. (56) 


(49),(60) and (61) give a closed form expression for the 
Laplace transform of the dispersion tensor at all 

times and is the key result of the paper. 


IV. THE EFFECTIVE DIFFUSION TENSOR 


Here we extract the late time effective diffusion tensor 


Dij. The Eqs. (60) and (611 are written in Laplace space 


and the large time limit is thus extracted from the small 
s behavior, for which 


CTy(s) ~ 


2 D, 


(s ^ 0). 


(62) 
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By definition, 

poo 

P'(x,s|y)=/ dt e"®*[P(x,t|y) - Po(x)]. (63) 

Jo 

Taking s = 0 in the above expression, we find that P' 
is finite at s = 0 (as P(x,t|y) —Po(x) as t —)• oo). 
Therefore, we can define a function G such that 

poo 

P(x,0|y) = / dt [P(x,t|y) -Po(x)] = G(x|y). (64) 
Jo 

Applying the operator H to the above equality, it is clear 
that G obeys 

PxG(x|y) =(5(x-y)-Po(x), (65) 

and also from the conservation of probability G satisfies 

[ dx G(x|y) = 0. (66) 

Jn 

The function G is therefore the pseudo-Green’s function 
[34] of the operator H (that is, the inverse of H in the 
subspace orthogonal to the uniform function). 

Since P'(x, s = 0|y) is finite, we see that the inverse 
Laplace transform at late times can be extracted from 
Eqs. (60) and (611 upon setting s = 0 in P'. This imme¬ 
diately yields 

Dij = / dx Kij(x)Po(x) 

Jn 

+ ^ ^ C^X [Mj(x)/j(x) -h Wj (x)/,(x)] 


'Sn 


dS'i(x)[K,i(x)/j(x) -h Kj;(x)/*(x)], (67) 


with 


/j W=- f -Po(y)G(x|y)M*(y) 
Jn 


+ / c^5'/(y)Po(y)G(x|y)/Cij(y). (68) 

JSa 

For computational purposes, it is useful to find the partial 
differential equation satisfied by fi. Acting with H on the 
above expression and using Eq. (65) leads to 


Px/*(x) =-M-(x)Po(x)-hPo(x) [ dy [M-(y)Po(y)] 

Jn 

--Po(x) [ dSi{y)Po{y)Kuiy), (69) 

JSn 

which is valid for x S (outside the obstacles). An 
equivalent equation for fi is 

Px/»(x) = -M-(x)Po(x) 

-fPo(x) f dy {u*(y)Po(y) - (9yjPo(y)Kfcj(y)]}. (70) 
Jn 


This equation for /i(x) (for x outside the obstacles) is 
supplemented by the following conditions. First, because 
of the conservation of probability, it follows from Eq. (66) 
that 


[ dx. /,(x) = 0. 
Jn 


(71) 


Second, /i(x) is periodic on fl. Third, at the surface 
obstacles, fi satisfies 

nj{uj{x)f,{x) - 9,,jKfcj(x)/,(x)]} = -njKj^{x)Po{x) 

(72) 

It is actually not obvious to derive the boundary condi¬ 
tions (|72|) from Eq. (1^. The validity of these boundary 


conditions is checked explicitly in Appendix]^ where we 
demonstrate that Eq. (|6^ is the actual solution of the 
Eqs. ( |70|71|72D , proving that our formulation is correct. 

One may alternatively express the diffusion tensor by 
defining a function f* as 


/*(y) = - [ dx Po(y)G(x|y)M,(x) 
Jn 


1 

2 Jn 
\ ^ 


in 

+ [ dSi{x)Po{y)G{x\y)Ku{x), (73) 

JSn 

in terms of which the effective diffusion tensor reads 

D^J = dy Kij(y)Po(y) 

Jn 

dy K(y)/j‘(y) + u*{y)f*{y)] 
dSi{y)[Ka{y)fj[y) + Hji{y)f*{y)]- (74) 

Obviously, 

/*(y)=- [ dx Po(x)G*(y|x)Mi(x) 

Jn 

P [ dSi{x)Po{x)G{y\x)Kii{x), (75) 

JSa 

where G* is the pseudo-Green function associated for the 
operator H* of the time reversed process, with the drift 
field u*, that reads 

77y(-) = 9yJ<(y)(-)] - dy,dy^[K^j{y){-)] (76) 

Hence, f* satisfies the same equations as fi if one replaces 
Ui by u* and vice-versa. As a consequence, /* is the 
solution of 

Hyf*(y) = - uziy)Poiy) + Po(y) [ dx [u,(x)Po(x)] 

Jn 

- Po{y) [ dSi{x)Poix)Kii{x) (77) 

JSn 

with periodic boundary conditions and also with the con¬ 
dition 

nj{ujiy)fHy) - dy^[ii-kj{y)fHy)]} = -noK,,i{y)Po{y) 

(78) 









at the surface of the obstacles, together with the orthog¬ 
onality condition 


/ dy/;(y) = 0. (79) 

jn 

Thus, we have two formulas for the computation of the 
effective diffusion tensor. 

Having two formulations for the Kubo formula is a 
useful check on numerical solutions as the two resulting 
results for the diffusion tensor can be compared. Fur¬ 
thermore, we will now see that both the fields fi and 
f* are needed to compute the leading order finite time 
correction to the time dependent diffusion tensor. 


V. LATE TIME CORRECTIONS TO THE 
EFFECTIVE DIFFUSION TENSOR 


In section IV we have extracted, from the full time 
dependent Kubo formula, the effective late time limit of 
the diffusion tensor Dij. While many examples of such 
formulas exist for special cases, little is known about how 
the time dependent diffusion tensor relaxes to its asymp¬ 
totic limit. To our knowledge the only examples known 
are for diffusion in a periodic potential m and diffusion 
in a system with periodic diffusivity [SS]. Both of these 
examples have steady states with zero current, here we 
will extend these results to the most general cases both 
with and without current. Here we extract the late time 
correction tensor Cij that describes the approach to the 
diffusive limit of the system. 


Expanding Eq. (63) in powers of s, we obtain 


P'(x, s|y) ~ G(x|y) + s G(2)(x|y) + ... (80) 


with 


pOO 

G(2)(x|y) = - dtt [P(x,t|y) - Po(x)]. 
JQ 


(81) 
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The correction tensor G^ is identified from the relation 
for small s 


+ s Cij + ...) 


(s^O). (85) 


Inserting this relation into Eqs. (60) and ( |61[ ), and using 
([sol) and ([84|) leads to 


G„ = - ^ 
2 


dz 


{ rfx[u,(x)/j(z) -F Mj(x)/,(z)]G(x|z) 

-[ dS'i(x)[K,i(x)/j(z)-f Kj7(x)/j(z)]G(x|z)|, (86) 

JSn ^ 

where we have used the definition of / in Eq. ( [6^ to per¬ 
form the integration over y. Using the definition of f* in 
Eq. (75) and the relation G*(y|x)Po(x) = G(x|y)Po(y): 
we finally obtain 


^ Po{^) 


This relation is a compact and explicit form for the cor¬ 
rection tensor Cij and is the main result of this section. 
One can also derive the above result by decomposing P' 
in terms of left and right eigenvector of the Fokker-Planck 
operator H. In doing this is is straightforward to see 
that the temporal corrections at the next order decay as 
exp(—Alt), where Ai is the lowest positive eigenvalue of 
H which will be strictly positive given that the domain 
n is taken to be finite. 


VI. SPECIAL CASES AND COMPARISON 
WITH EXISTING RESULTS 

Here we examine the form the Kubo formula takes for a 
number of systems and compare them to existing results 
in the literature. 


A. Flow in frozen incompressible velocity fields 
with isotropic constant diffusion tensor 


Since the motion of the tracer particles is Markovian 
(memoryless), we can write the equality 

-P(x,t|y)=/ dz P(x,t-fi|z)P(z,ti|y), (82) 
Jn 


The Fokker-Planck transport operator for a tracer ad- 
vected by an incompressible velocity field v with constant 
isotropic diffusivity Kii(x) = KoSij is defined via 

-ffx = -Koda;idxi + Vi{-K)dx^ . (88) 


which holds for 0 < ti < t. Integrating this relation over 
ti leads to 

tP(x,t|y)= f dti /" dz P(x,t - ti|z)P(z,ti|y). 

Jo dn 

(83) 

Substracting tPg (x) and taking the Laplace transform for 
small values of s, we obtain 

G(2)(x|y) = - /'dzG(x|z)G(z|y). (84) 

Jn 


Here the drift field is thus Ui = Vi. From the incom¬ 
pressibility condition dx^Vi = 0, we see that the steady 
state distribution on the unit cell H is uniform and thus 
given by Po(x) = 1/|H|, where |H| denotes the available 
volume of the unit cell. The steady state current is thus 
given by Ji(x) = Ui(x)/|H|. We also find that u* = —Ui. 
Applying Eq. (67) leads to 


Ai = + l f d.x [fi{x)vj{x) + fj{x)v,{x)] 

^ Jn 

~Y {is W + c(A(x)/i(x)|. (89) 
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The equation satisfied by fi is a simplified form of Eq. 


The boundary conditions on the surface of obstacles are 
[see Eq. @] 


rijdxji = ni/\fl\, 


(91) 


(where we have assumed that no fluid enters the ob¬ 
stacles, ViTii = 0). It is then obvious that the last 
four expressions recover those in the hydrodynamics- 
homogenization literature (see e.g. Refs. [lll[35]), if we 
identify the function fi of this paper with —/i/jnj to 
match with the notations of Ref. [3S] . 


B. Systems with no steady state current 

A large class of models studied in statistical mechan¬ 
ics, such as diffusion in a periodic potential or diffusion in 
a medium having periodic diffusivity, have no current in 
the steady state. In these cases the original and its time- 
reversed processes have the same Eokker-Planck evolu¬ 
tion operator, and consequently the same pseudo-Green’s 
function. However in general we have that 

G(x|y)Po(y) = Ho(x)G*(y|x), (92) 

where G*(y|x) denotes the pseudo-Green’s function for 
H*. In the case where there is no current we know how¬ 
ever that G*(y|x) = G(x|y) and thus we find that the 
operator 


K{x,y) = G(x|y)Po(y) = Po(x)G(y|x) (93) 


is symmetric. In the currentless case we also have 


the diffusing particle becomes trapped in local minima of 
the potential. This slowing down due to drift also occurs 
for diffusion in a medium of spatially varying diffusivity 
where 


(96) 

here Pq = 1/|H| and Ji = 0. The drift field in this case 
is, in our notation, given by 


Wi(x) = 92,^-Kij(x), (97) 

and in terms of the pseudo-Greens’ function G we find 

- [dx^Kik{x)][dx,Ku{y)]G{x\y), (98) 


which recovers the formula derived in Ref. [H]. We see 
that the first term in the above is the spatial average of 
the local diffusivity tensor while the second terms gives 
a negative contribution when i = j. Einally we remark 
that the case of diffusion with constant diffusivity kq in 
a periodic potential $ gives Ui = — ko/35x,$(x), where f3 
is the inverse temperature. Here the steady state distri¬ 
bution Pq is given by the Gibbs-Boltzmann distribution 
on the unit cell H 


exp [-/3$(x)] 
G)ixj -- — -, 


(99) 


where Zq is the partition function normalizing the dis¬ 
tribution, thus recovering the result of Ref. m for this 
particular case. 


VII. GENERAL RESULTS IN ONE DIMENSION 


Dij = / dx Kij(x)Po(x) 

Jq 

- ^ dxdy[fo(x)Mj(y) -h u^(x)fo(y)]A:(x, y), (94) 


from Eq. (10). The operator K, being symmetric, must 


have real eigenvalues, however K has the same eigenval¬ 
ues as G which must have eigenvalues with a positive 
real part if a steady state regime can be attained. Con¬ 
sequently under these assumptions K must be a positive 
operator and we must have that the change in the diag¬ 
onal terms of the diffusion tensor due to the presence of 
drift or inclusions is negative, i.e. 


Pii (^iz)o — 



dxdyui{x.)ui{y)K{x,y) < 0, (95) 


that is to say that diffusion in equilibrium systems is 
slowed down by the presence of the drift. This is physi¬ 
cally obvious for say diffusion in periodic potentials where 


We consider a generic one dimensional system with 
periodicity denoted by L. Here we can give explicit for¬ 
mula for all terms in the static Kubo formula, although 
the computation is surprisingly involved. In what follows 
we will use notation similar to Refs. [551 HH] to facili¬ 
tate comparison with their results. In one dimension the 
steady state current is constant in space and the steady 
state probability distribution is given by 

Po{x) = JI+{x), (100) 

where the term 1+ reads 

I+ix) = ^^P(^(^(( f dx' exp{-T{x')), (101) 

k[x) Jx 

with 

r{x)= Tdx'^^. ( 102 ) 

Jo 
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Due to the periodicity of u and k the function F obeys 
the relation 


r(a; + L) =r(x) +r(L). 


(103) 


When r(L) = 0 the system clearly has a steady state 
equilibrium distribution with no current. In writing Eq. 
(Toll) we have assumed, without loss of generality, that 
r(I)> 0 so that the integral on the right hand side con¬ 
verges. The steady state current is then obtained from 
the condition of normalization of Pq and is thus given by 


J = 


fo dx I+{x) 


(104) 


and consequently the effective drift is given hy V = JL. 

To compute the effective diffusion constant we use the 
representation for D given in Eq. (741 in its one dimen¬ 
sional version, that is to say 

pL pL 

D= dx k{x)Po{x) + / dxf*{x)u*{x). 

Jo Jo 

We now define f*{x) = g{x)Po{x) and write the diffusion 
constant as 

D = [ dx k{x)Po{x)+ [ dxg{x) IJ —— [K(a;)Po(2;)] 
Jo Jo I dx 

Defining g(x) = —x + h{x) we find that h obeys 


dJh dh 
K-r^ + u— = JL. 
dx^ dx 


A first integration of this equation for h yields 
dh{x) 


dx 


= JLI- (x), 


where 


/-(a;) = exp(-r(a;)) / 

J —c 


dx' 


, ex p(r(a:')) 

k{x') 


(105) 


(106) 


(107) 


and where we have used the fact that dg/dx = — 1+dh/dx 
is periodic and thus dh/dx must be periodic. Integrating 
again yields 


h{x) = JoL[ dx'I_{x') + C 

Jo 


(108) 


where C is an integration constant that will 
be determined from the orthogonality condition 
dx Po( x)g( cc) = 0. It is not immediately obvious 


from Eq. ( 108| ) that g{x) = —x + h is periodic, however 
it is easy to see that 


g{x + L) = g{x) — L + JL f dx I-{x). 

Jo 

Now by integration by parts one finds 

pL pL 

/ dx I-{x) = dx /+(a:), 

Jo Jo 


(109) 


( 110 ) 


and using this and Eq. (1041 in Eq. (109) yields g{x) = 


g{x + L). Incidentally this shows that we may also write 

1 


J = 


fo dx I-(x) 


( 111 ) 


To proceed it is useful to define the two new functions 


A+(x) = 

r dx' 

7-00 «(P) ’ 

(112) 

A_(x) = 

f dx'exp (—r(x')). 

J X 

(113) 

Eurther more one can 

show that 


A_|_ (x H- Lj) 

= exp(r(L)) A+(x), 

(114) 

A- (x -I- L) 

= exp (—r(L)) A_(x). 

(115) 


In terms of these two functions we can write /+ and /_ 
as 


^ , , dA-{x) . , , 

/_(x) =- ^A+{x). 


(116) 

(117) 


Computing the constant C in Eq. (108) using the orthog¬ 
onality relation then yields 


L>= dx 

Jo 


k{x)Po{x)^ 


-I- J [—X -|- h{x) + LxPo{x) — Lh{x)Po{x)\ 


(118) 


where we have carried out an integration by parts in the 
last term of Eq. (105). The function h is explicitly given 
by 


h{x) = JL f dx' I-(x'), 

Jo 


(119) 


and using the expression Eq. (117) for /_ it can be writ¬ 
ten as 


h{x) = — JL 


A-|_(x)A_(x) — A+(0)A_(0) 

./ Mx') 


- / dx' 

Jo 


J 


( 120 ) 


We also, by integration by parts, note the identity 

pL px pL 

dx dx' Po{x') = L — / dx xPo{x). (121) 

Jo Jo Jo 


Using these relations in Eq. (118), after some algebra, 
we obtain the compact expression 

L'^ ffj dx K{x)I±{xYlzn{x) 

D= ^ V / ±w j,^22) 


Jo dx /±(x)3 
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where ± indicates that the index can be taken to be + or 
—. One can verify that Eq. ( 122\ agrees with the results 
of Reimann et al. [5S1 [55] for the case of diffusion in a 
tilted potential (j){x) = V{x) — fx, where V is periodic. 
It also agrees with the formulas presented in Ref. m 
where the local mobility is also varying. The result given 
here shows how the diffusion constant of any one dimen¬ 
sional system with periodic diffusivity and drift can be 
obtained. 


VIII. CONCLUSIONS 

The problem of computing effective transport coeffi¬ 
cients occurs in many settings of fundamental physics and 
has many applications. Many results exist in the litera¬ 
ture, for example in fluid mechanics using homogeniza¬ 
tion theory pT)l - l20] as well as in statistical physics [2TH25] 
where both exact and approximate results exist for prob¬ 
lems such as diffusion in periodic or random potentials 
and diffusion in media with locally varying diffusivity. 
In these latter studies most results have been obtained 
for systems with equilibrium (currentless) steady states, 
given for example by the Gibbs Boltzmann distribution 
for diffusion in a periodic potential. 

More recently results have been derived for diffusion in 
systems with non zero currents [26II32j In such non equi¬ 
librium systems interesting new physics such as massively 
increased diffusivity due to an external applied force has 
been discovered. This enhancement occurs, in the weak 
noise limit, at the external field where all minima disap¬ 
pear from the overall potential to be replaced by points 
of inflection |5^. The difference between different tra¬ 
jectories at this critical point is enhanced and dispersion 
in thus increased. Using the Kubo formula, given here, 
one could investigate the effect of an externally applied 
field on diffusion in higher dimensional periodic poten¬ 
tials. The formula given here cannot in general be eval¬ 
uated analytically in higher than one dimension, how¬ 
ever their numerical evaluation is straightforward and 
the transport coefficients can be accurately determined 
using standard packages to solve partial differential equa¬ 
tions. Of course many of the systems considered here 
can be studied via numerical simulations where one in¬ 
tegrates the corresponding SDE, however the simulation 
approach suffers from the need to estimate the errors due 
to statistical fluctuations and in the case of obstacles the 
correct imposition of the no-flux boundary condition re¬ 
quires careful treatment and is far from being obvious 
[40H42|. 

Many other problems can be tackled using the present 
approach, in Ref. [33j it was shown how the presence 
of an externally applied uniform field on a system with 
spatially varying diffusivity modifies the dispersion of a 
cloud of tracer particles. It was shown that the diffu¬ 
sion in the direction of the applied force could be hugely 
increased for systems in two or more dimensions. In ad¬ 
dition it was shown that the components of the effective 


diffusion tensor can exhibit counter intuitive monotonic 
behavior. The results thus imply that one can control 
the dispersion of a diffusing cloud with an external field, 
and such an effect may well have useful applications. 

The formula given here are extremely general and 
should be valuable for the interpretation of experiments. 
For instance, potential fields in which colloidal particles 
diffuse are often generated by laser [43ll46] , in which case 
absorption of laser light can also lead to temperature gra¬ 
dients. The variation of the local temperature will clearly 
influence local transport properties through the tempera¬ 
ture dependence of the surrounding fluid viscosity as well 
as due to thermodynamic forces associated with temper¬ 
ature gradients, i.e the Soret effect 

Finally there has been much recent study of fluctuation 
dissipation relations (FDR) in non-equilibrium systems 
[37HS0I- A notable example of this type of FDR is the 
Stokes-Einstein relation between the differential response 
of the mean velocity Vi and the diffusion tensor. In our 
preliminary study [53], we showed how the presence of 
a steady state current leads to corrections to the usual 
Stokes-Einstein formula. However the general analysis 
of the correction term needs further study in order to 
fully understand the non-equilibrium physics it encodes, 
as well as to make contact with the literature on the sub¬ 
ject [571150] . In the literature in question results are most 
often given in terms of an integral over a time dependent 
violation factor, the results given here could be useful 
to understand these results in terms of the static, time 
independent quantities, used here. 


Appendix A: Details on the numerical simulations 


Here we briefly describe the algorithm used to simulate 
the motion of a Brownian walker in a 2D array of circular 
disks, in the presence of a force F, leading to the results 
presented in Figj^b). Let X be the position at time t 
of the walker, then we define the unit vector with 
direction X — Xc, with Xc the nearest disk center, 
makes an angle a with the x—direction. During the time 
step At, the motion in the direction is modified by 
the increment 

dX± = '/Kifi{r±/VAt) +uiif2ir±/VAi) + h cos a At 

(Al) 


Here h = PFD, r± is the distance to the nearest disk 
surface, and /i and /2 are the functions calculated by 
Peters et al. in the case of reflecting boundaries m, and 
U|| = ±^DAt with probabilities 1/2. In the direction 
parallel to the nearest surface obstacles, the increment is 
dXu = —hsinaAt u, where u is a stochastic variable 


with zero mean and variance (u^) = 2D At. With this 
algorithm, when the particle is close to the boundaries, 
r_L ^ VAt, the terms in fi and /2 dominate the incre¬ 
ment in the perpendicular direction in Eq. (Al) if At 
is small enough. These terms account for the fact that. 
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since trajectories that cross the surface of the obstacles 
are forbidden, the averaging over authorized trajectories 
between t and t + At gives rise to an effective drift ori¬ 
ented towards the exterior of the obstacles, described by 
the term /i, and also modifies the variance of these tra¬ 
jectories (term / 2 ), see Ref [40] for details. If the particle 
is far from the obstacles, rj_ 3> then /i —>■ 0 and 

/2 —t and one recovers Brownian motion under force 
in the absence of obstacles. Although the algorithm of 
Ref |20| was proposed only for one-dimensional geome¬ 
tries, we expect that it is also valid in the present 2D 
case, since for sufficiently small At the curvature of the 
surface should play no role. For the results of Fig. Bb). 
the time step was At = 10 /D and was checked to 
be sufficiently small so that convergence was reached for 
the calculation of which were estimated over 500,000 
trajectories. The results of simulations match with the 
predictions of our (exact) Kubo formulas. 


The following condition also holds 


dy Po(y)G(x|y) = 0. 


(B4) 


Using the above relations and Eqs. (69) and 0 we de¬ 
duce that 


Mi{x) = fi{x) + [ dy Po(y)G(x|y)u-(y). (B5) 

Jn 


On the other hand, using the explicit expressions of H 
[Eq.([^] and integrating by parts two times the term con¬ 
taining H in Eq. (Bl), we obtain 


Appendix B: Boundary condition (72| for the 
function fi 

In this appendix, we show that the formulation of the 
problem under the form of Eqs. ( |70|71|72D is correct by 
checking that Eq. (68) is the actual solution of these equa¬ 


tions. We consider the quantity 


M,(x)= / dy {[ij1;G(x|y)]/i(y) - G(x|y)[i7y/,(y)]} 
Jn 

(Bl) 


where the adjoint operator reads 

Hi = -Ui{y)dy, - K^j{y)dyMy^. 


(B2) 


The pseudo-Green function G can be shown to satisfy 
the adjoint equation 


id^G(x|y) = d(x - y) - Po(x), 


(B3) 


Mi(x) = 

dSj{y)G{x\y){-Uj{y)fiiy) + dy^ K/c(y)/z(y)]} 


ISu 


>Sn 


dSjiy)fiy)Kkjiy)dy^G{x\y). 


(B6) 


The last term of this equation vanishes, since the relation 
dSjKkj{y)P{x,t\y) = 0 is known for reflecting bound¬ 
aries [3|. Then, using the boundary conditions (72), we 
obtain 


M,(x) = f dSj{y)G{x\y)Kjkix), (B7) 
JSn 


which, when compared to Eq. (B5), is exactly Eq. ( 68), 
thereby proving that the formulation of Eqs. ( |70i72|7lD 
is correct. 
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